home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / fft / real_main.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-08-22  |  3.6 KB  |  146 lines
  1. /* fft/real_main.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <math.h>
  23.  
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_complex.h>
  26. #include <gsl/gsl_fft_real.h>
  27.  
  28. #include "real_pass.h"
  29.  
  30. int
  31. FUNCTION(gsl_fft_real,transform) (BASE data[], const size_t stride, const size_t n,
  32.                   const TYPE(gsl_fft_real_wavetable) * wavetable,
  33.                                   TYPE(gsl_fft_real_workspace) * work)
  34. {
  35.   const size_t nf = wavetable->nf;
  36.  
  37.   size_t i;
  38.  
  39.   size_t q, product = 1;
  40.   size_t tskip;
  41.   size_t product_1;
  42.  
  43.   BASE *const scratch = work->scratch;
  44.   TYPE(gsl_complex) *twiddle1, *twiddle2, *twiddle3, *twiddle4;
  45.  
  46.   size_t state = 0;
  47.   BASE *in = data;
  48.   size_t istride = stride ;
  49.   BASE *out = scratch;
  50.   size_t ostride = 1 ;
  51.   
  52.   if (n == 0)
  53.     {
  54.       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
  55.     }
  56.  
  57.   if (n == 1)
  58.     {                /* FFT of one data point is the identity */
  59.       return 0;
  60.     }
  61.  
  62.   if (n != wavetable->n)
  63.     {
  64.       GSL_ERROR ("wavetable does not match length of data", GSL_EINVAL);
  65.     }
  66.  
  67.   if (n != work->n)
  68.     {
  69.       GSL_ERROR ("workspace does not match length of data", GSL_EINVAL);
  70.     }
  71.  
  72.   for (i = 0; i < nf; i++)
  73.     {
  74.       const size_t factor = wavetable->factor[i];
  75.       product_1 = product;
  76.       product *= factor;
  77.       q = n / product;
  78.  
  79.       tskip = (product_1 + 1) / 2 - 1;
  80.  
  81.       if (state == 0)
  82.     {
  83.       in = data;
  84.       istride = stride;
  85.       out = scratch;
  86.       ostride = 1;
  87.       state = 1;
  88.     }
  89.       else
  90.     {
  91.       in = scratch;
  92.       istride = 1;
  93.       out = data;
  94.       ostride = stride;
  95.       state = 0;
  96.     }
  97.  
  98.       if (factor == 2)
  99.     {
  100.       twiddle1 = wavetable->twiddle[i];
  101.       FUNCTION(fft_real,pass_2) (in, istride, out, ostride, product, n, twiddle1);
  102.     }
  103.       else if (factor == 3)
  104.     {
  105.       twiddle1 = wavetable->twiddle[i];
  106.       twiddle2 = twiddle1 + tskip;
  107.       FUNCTION(fft_real,pass_3) (in, istride, out, ostride, product, n, twiddle1,
  108.                    twiddle2);
  109.     }
  110.       else if (factor == 4)
  111.     {
  112.       twiddle1 = wavetable->twiddle[i];
  113.       twiddle2 = twiddle1 + tskip;
  114.       twiddle3 = twiddle2 + tskip;
  115.       FUNCTION(fft_real,pass_4) (in, istride, out, ostride, product, n, twiddle1,
  116.                      twiddle2, twiddle3);
  117.     }
  118.       else if (factor == 5)
  119.     {
  120.       twiddle1 = wavetable->twiddle[i];
  121.       twiddle2 = twiddle1 + tskip;
  122.       twiddle3 = twiddle2 + tskip;
  123.       twiddle4 = twiddle3 + tskip;
  124.       FUNCTION(fft_real,pass_5) (in, istride, out, ostride, product, n, twiddle1,
  125.                      twiddle2, twiddle3, twiddle4);
  126.     }
  127.       else
  128.     {
  129.       twiddle1 = wavetable->twiddle[i];
  130.       FUNCTION(fft_real,pass_n) (in, istride, out, ostride, factor, product, n,
  131.                      twiddle1);
  132.     }
  133.     }
  134.  
  135.   if (state == 1)        /* copy results back from scratch to data */
  136.     {
  137.       for (i = 0; i < n; i++)
  138.     {
  139.       data[stride*i] = scratch[i] ;
  140.     }
  141.     }
  142.  
  143.   return 0;
  144.  
  145. }
  146.